All runs Hypnea musciformis PHOTOSYNTHESIS Analysis, Script Chunks, and Plots

This is the analysis of all the Hypnea musciformis salinity and nutrient experiments conducted on the lanai in St. John 616 from September 2021 to October 2022. These experiments incorporated four paired salinity and nutrient treatments with three temperatures. Each of the first four runs produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled in February, April, and October 2022. This output reflects all data totaling six treatments for Hypnea only.

Packages loaded:

library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
#library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load and prepare the dataset

Open the output dataset generated by the ps_script_clean_to_ek_alpha.R script in the phytotools_alpha_ek project This file was normalized to quantum efficiency of photosynthesis as this seems to be more accurate (changed in fitWebb input) per Silsbe and Kromkamp (2012)

all_runs_photosyn_data <- read.csv("../data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")

Assign run as a factor

all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)

Assign temperature as a factor

all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)

Assign treatment as characters from integers then to factors

all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))

Assign deltaNPQ as a factor

all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)

Subset the data and toggle between the species for output. Use Day 9 for final analysis ONLY This will also assign the proper labels for plots There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21

hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

Add a column for growth rate from growth rate dataset to the already subsetted hypnea data frame

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100

Subset for hypnea only and calculate growth rate from final and initial weights Removed one dead individual hm6-4 from 10/9/21 (hm6-4) and another dead from 11/12/21 (also (hm6-4)

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

#Run the model Run model for rETRmax with two fixed effect variables and three random effects variables

#run model without interaction between the treatments and temperature
all_runs_photosyn_model_hyp <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run)
                                      + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)

rETRmax – Make a histogram and residual plots of the data

hist(hypnea$rETRmax, main = paste("Hypnea musciformis rETRmax"), col = "maroon", labels = TRUE)

#or
hypnea %>% ggplot(aes(rETRmax)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

plot(resid(all_runs_photosyn_model_hyp) ~ fitted(all_runs_photosyn_model_hyp))

qqnorm(resid(all_runs_photosyn_model_hyp))
qqline(resid(all_runs_photosyn_model_hyp))

rETRmax – Check the performance of the model

performance::check_model(all_runs_photosyn_model_hyp)

These outputs show the model is acceptable

rETRmax – Check r2 for model fit and print the model statistics summary

r.squaredGLMM(all_runs_photosyn_model_hyp)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##            R2m       R2c
## [1,] 0.3005666 0.6509689

summary(all_runs_photosyn_model_hyp)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: hypnea
## 
##      AIC      BIC   logLik deviance df.resid 
##   2334.8   2378.7  -1155.4   2310.8      274 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2015 -0.5040 -0.0661  0.4101  4.0765 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  99.958   9.998  
##  Run       (Intercept)  16.273   4.034  
##  RLC.Order (Intercept)   1.092   1.045  
##  Residual              116.863  10.810  
## Number of obs: 286, groups:  Plant.ID, 143; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    67.4755     3.9362   8.3868  17.142 8.02e-08 ***
## Treatment1      0.2121     4.5415   8.1999   0.047    0.964    
## Treatment2     -1.7598     4.5415   8.1999  -0.387    0.708    
## Treatment2.5    9.9921     5.7830   4.6952   1.728    0.148    
## Treatment3      3.9787     4.5415   8.1999   0.876    0.406    
## Treatment4      1.9950     4.5554   8.3036   0.438    0.673    
## Temperature27 -19.1213     2.8560  11.6971  -6.695 2.51e-05 ***
## Temperature30 -20.0265     2.8042  21.9622  -7.142 3.72e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.709                                          
## Treatment2  -0.709  0.882                                   
## Treatmnt2.5 -0.556  0.482  0.482                            
## Treatment3  -0.709  0.882  0.882  0.482                     
## Treatment4  -0.708  0.879  0.879  0.480  0.879              
## Temperatr27 -0.359  0.005  0.005  0.000  0.005  0.005       
## Temperatr30 -0.357  0.000  0.000  0.000  0.000  0.004  0.494

Use Bartlett’s test to check for equal variance

bartlett.test(rETRmax ~ Treatment, data = hypnea)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  rETRmax by Treatment
## Bartlett's K-squared = 16.492, df = 5, p-value = 0.005571

Run Welch’s ANOVA if not equal variances

welch_anova_treatment <- oneway.test(rETRmax ~ Treatment, data = hypnea, var.equal = FALSE)
welch_anova_treatment

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  rETRmax and Treatment
## F = 2.307, num df = 5.00, denom df = 130.17, p-value = 0.04796

welch_anova_temp <- oneway.test(rETRmax ~ Temperature, data = hypnea, var.equal = FALSE)
welch_anova_temp

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  rETRmax and Temperature
## F = 35.12, num df = 2.00, denom df = 183.65, p-value = 1.214e-13

games_howell_test(hypnea, rETRmax ~ Treatment, conf.level = 0.95, detailed = FALSE)

rETRmax – Plot and make a table of the results. Also get means for the treatments

plot(allEffects(all_runs_photosyn_model_hyp))

tab_model(all_runs_photosyn_model_hyp)

	rETRmax
Predictors	Estimates	CI	p
(Intercept)	67.48	59.73 – 75.22	<0.001
Treatment [1]	0.21	-8.73 – 9.15	0.963
Treatment [2]	-1.76	-10.70 – 7.18	0.699
Treatment [2.5]	9.99	-1.39 – 21.38	0.085
Treatment [3]	3.98	-4.96 – 12.92	0.382
Treatment [4]	2.00	-6.97 – 10.96	0.662
Temperature [27]	-19.12	-24.74 – -13.50	<0.001
Temperature [30]	-20.03	-25.55 – -14.51	<0.001
Random Effects
σ²	116.86
τ₀₀ _Plant.ID	99.96
τ₀₀ _Run	16.27
τ₀₀ _RLC.Order	1.09
ICC	0.50
N _Run	8
N _Plant.ID	143
N _RLC.Order	6
Observations	286
Marginal R² / Conditional R²	0.301 / 0.651

hypnea %>% group_by(Treatment) %>% summarise_at(vars(rETRmax), list(mean = mean))

hypnea %>% ggplot(aes(treatment_graph, rETRmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "rETRmax (μmols electrons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +   ylim(-1, 170) + stat_mean() + 
  geom_hline(yintercept=0, color = "purple", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

Plot a regression between the photosynthetic independent variables of interest and growth rate

#rETRmax vs. Growth rate
hypnea_growth_etr_graph <- ggplot(hypnea, aes(x=rETRmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
hypnea_growth_etr_graph

## `geom_smooth()` using formula 'y ~ x'

#rETRmax vs. Ek
hypnea_etr_ek_graph <- ggplot(hypnea, aes(x=rETRmax, y=ek.1)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis rETRmax vs Ek", x = "rETRmax (μmols electrons m-2 s-1)", 
       y = "Ek (μmols photons m-2 s-1)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
hypnea_etr_ek_graph

## `geom_smooth()` using formula 'y ~ x'

Temperature did not have a significant effect on the outcome for rETRmax

Ek – Run the model

Run model for minimum saturating irradiance (Ek) with two fixed effect variables and three random effects variables

all_runs_photosyn_model_ek <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run)
                                      + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea)

Ek – Make a histogram and residual plots of the data for hypnea

hist(hypnea$ek.1, main = paste("Hypnea musciformis Ek"), col = "darkolivegreen3", labels = TRUE)

plot(resid(all_runs_photosyn_model_ek) ~ fitted(all_runs_photosyn_model_ek))

qqnorm(resid(all_runs_photosyn_model_ek))
qqline(resid(all_runs_photosyn_model_ek))

hypnea %>% ggplot(aes(ek.1)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

Ek – Check the performance of the model

performance::check_model(all_runs_photosyn_model_ek)

Ek – Check r2 for model fit and print the model statistics summary

r.squaredGLMM(all_runs_photosyn_model_ek)

##            R2m       R2c
## [1,] 0.3740684 0.5666269

summary(all_runs_photosyn_model_ek)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: hypnea
## 
## REML criterion at convergence: 2624.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0840 -0.5727 -0.0771  0.4386  3.8792 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept) 101.41   10.07   
##  Run       (Intercept) 127.93   11.31   
##  RLC.Order (Intercept)  13.61    3.69   
##  Residual              546.79   23.38   
## Number of obs: 286, groups:  Plant.ID, 143; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     76.746      9.327   6.547   8.229  0.00011 ***
## Treatment1      18.127     10.821   6.094   1.675  0.14415    
## Treatment2      17.454     10.821   6.094   1.613  0.15710    
## Treatment2.5    42.029     14.796   4.649   2.841  0.03936 *  
## Treatment3      17.248     10.821   6.094   1.594  0.16129    
## Treatment4      24.509     10.845   6.150   2.260  0.06348 .  
## Temperature27  -35.835      4.845  20.376  -7.397 3.41e-07 ***
## Temperature30  -39.110      4.540  48.669  -8.615 2.35e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.768                                          
## Treatment2  -0.768  0.903                                   
## Treatmnt2.5 -0.561  0.484  0.484                            
## Treatment3  -0.768  0.903  0.903  0.484                     
## Treatment4  -0.767  0.901  0.901  0.483  0.901              
## Temperatr27 -0.251  0.002  0.002  0.000  0.002  0.000       
## Temperatr30 -0.245  0.000  0.000  0.000  0.000  0.005  0.480

Ek – Run Bartlett’s test and Welch’s ANOVA with Games Howell test for pairwise comparisons

bartlett.test(ek.1 ~ Treatment, data = hypnea)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  ek.1 by Treatment
## Bartlett's K-squared = 6.2901, df = 5, p-value = 0.279

welch_anova_treatment <- oneway.test(ek.1 ~ Treatment, data = hypnea, var.equal = FALSE)
welch_anova_treatment

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  ek.1 and Treatment
## F = 9.3597, num df = 5.00, denom df = 130.34, p-value = 1.221e-07

welch_anova_temp <- oneway.test(ek.1 ~ Temperature, data = hypnea, var.equal = FALSE)
welch_anova_temp

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  ek.1 and Temperature
## F = 44.989, num df = 2.00, denom df = 187.18, p-value < 2.2e-16

games_howell_test(hypnea, ek.1 ~ Treatment, conf.level = 0.95, detailed = TRUE)

Ek – Plots and tables for the results

tab_model(all_runs_photosyn_model_ek)

	ek.1
Predictors	Estimates	CI	p
(Intercept)	76.75	58.39 – 95.11	<0.001
Treatment [1]	18.13	-3.18 – 39.43	0.095
Treatment [2]	17.45	-3.85 – 38.76	0.108
Treatment [2.5]	42.03	12.90 – 71.16	0.005
Treatment [3]	17.25	-4.05 – 38.55	0.112
Treatment [4]	24.51	3.16 – 45.86	0.025
Temperature [27]	-35.84	-45.37 – -26.30	<0.001
Temperature [30]	-39.11	-48.05 – -30.17	<0.001
Random Effects
σ²	546.79
τ₀₀ _Plant.ID	101.41
τ₀₀ _Run	127.93
τ₀₀ _RLC.Order	13.61
ICC	0.31
N _Run	8
N _Plant.ID	143
N _RLC.Order	6
Observations	286
Marginal R² / Conditional R²	0.374 / 0.567

plot(allEffects(all_runs_photosyn_model_ek))

hypnea %>% group_by(Treatment) %>% summarise_at(vars(ek.1), list(mean = mean))

hypnea %>% ggplot(aes(treatment_graph, ek.1)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Ek (μmols photons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +   ylim(-1, 240) + stat_mean() + 
  geom_hline(yintercept=0, color = "purple", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

plot linear regression Ek with growth rate

hypnea_growth_etr_graph <- ggplot(hypnea, aes(x=ek.1, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
hypnea_growth_etr_graph

## `geom_smooth()` using formula 'y ~ x'

Mixed Model for Photosynthesis in Hypnea musciformis

Angela Richards Donà

10/31/2022

Load and prepare the dataset

Ek – Run the model